#include <algorithm/analytical_model/algorithm/cuda/GeometryUtils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/SimplifiedMethod.cuh>
#include <algorithm/analytical_model/algorithm/cuda/Utils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/Utils.h>
#include <algorithm/analytical_model/algorithm/cuda/VectorMath.cuh>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/Data.h>
#include <algorithm/analytical_model/algorithm/precomp/NormalDistributionArray.h>
#include <utils/Format.h>

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <tuple>

#include <cuda_device_runtime_api.h>
#include <cuda_runtime_api.h>

using namespace solar::cuda;

namespace
{
    __device__ auto gnomonicProjection(float3 pixel_centre, float3 focus_centre, float3 helio_pos,
                                       float3 d_x, float3 d_y) -> float
    {
        // 计算接收器表面像素的四个顶点投影在成像平面上的坐标
        float3 v[4];
        float3 p_v[4];
        // 四个像素点
        v[0] = pixel_centre - d_x * 0.5 - d_y * 0.5;
        v[1] = pixel_centre + d_x * 0.5 - d_y * 0.5;
        v[2] = pixel_centre + d_x * 0.5 + d_y * 0.5;
        v[3] = pixel_centre - d_x * 0.5 + d_y * 0.5;
        // 投影在成像平面上

#pragma unroll
        for (int i = 0; i < 4; i++)
        {
            std::tie(p_v[i], std::ignore) =
                intersect(normalize(focus_centre - helio_pos), focus_centre, v[i],
                          normalize(v[i] - helio_pos));
        }
        // 计算面积
        float s = length(cross(p_v[0] - p_v[1], p_v[2] - p_v[1])) +
                  length(cross(p_v[0] - p_v[3], p_v[2] - p_v[3]));
        return 0.5F * s * rlength(cross(d_x, d_y));
    }

    __device__ auto projectToRecv(float3 pixel_center, float3 helio_center, float3 img_norm,
                                  float3 pixel_norm) -> float
    {
        float3 ray_dir = normalize(pixel_center - helio_center);
        return abs(dot(ray_dir, pixel_norm) / dot(ray_dir, img_norm));
    }

    __host__ __device__ constexpr auto computeAlpha(float u, int times) -> float
    {
        // float tmp = 0;
        // return std::modf(u * times, std::addressof(tmp));
        return u * times - int(u * times);
    }

    /**
     * @brief 每个线程按照id取接收器离散网格中心坐标并利用斜平行投影计算对应能量密度
     * @param d_array_ 记录光斑能量密度结果的数组
     * @param flux_param
     * 能量密度的系数，定日镜反射的太阳光总能量*反射光线与接收器的夹角余弦/2pi/sigma1/sigma2
     * @param row_col 接收器长宽对应的像素个数
     * @param focus_center 成像平面中心
     * @param center_bias 定日镜光斑在成像平面上的重心偏移，世界坐标系下
     * @param recv_vertex_0 接收器光斑图起始网格中心的世界坐标
     * @param d_x, d_y 为了获得网格点中心的世界坐标, recv_vertex_0 + i * d_x + j * d_y
     * @param x_axis 成像平面上的x轴
     * @param y_axis 成像平面上的y轴
     * @param z_axis 成像平面上的z轴
     */
    __global__ void bHFLCALMapFlux(float* d_array_, float flux_param, float sigma1, float sigma2,
                                   float correlation, int2 row_col, float3 focus_center,
                                   float3 center_bias, float3 recv_vertex_0, float3 d_x, float3 d_y,
                                   float3 x_axis, float3 y_axis, float3 z_axis, float3 plane_normal,
                                   float3 helio_pos)
    {
        long long my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= row_col.x * row_col.y)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        int2 pos_ij = make_int2(my_id % row_col.x, my_id / row_col.x);
        // 计算网格点的世界坐标
        float3 pos = recv_vertex_0 + pos_ij.x * d_y + pos_ij.y * d_x;
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        //    global_func::calcIntersection(z_axis, focus_center, pos, -z_axis, projected_pos);
        // 中心投影找对应点
        std::tie(projected_pos, std::ignore) =
            intersect(z_axis, focus_center, pos, normalize(helio_pos - pos));
        // 重心偏移
        //    projected_pos -= center_bias;
        //    global_func::calcIntersection(z_axis, focus_center, center_bias,
        //    normalize(helio_pos-pos), focus_center);
        float3 center_bias_image;
        std::tie(center_bias_image, std::ignore) =
            intersect(z_axis, focus_center, center_bias, -z_axis);

        // 投影点在成像平面上的二维坐标
        float u = dot(x_axis, projected_pos - center_bias_image);
        float v = dot(y_axis, projected_pos - center_bias_image);
        // 根据椭圆高斯分布计算格点处的能量
        // energy = E * f(x,y)
        float energy =
            flux_param * std::exp(-0.5 / (1 - correlation * correlation) *
                                  ((u * u) / (sigma1 * sigma1) + (v * v) / (sigma2 * sigma2) -
                                   2 * correlation * u * v / (sigma2 * sigma1)));
        // 乘投影系数
        // 斜平行投影Jacobian行列式
        //    energy *= parallelProjectionParam(x_axis,y_axis,z_axis,pos,helio_pos,plane_normal);
        //    energy *= abs(dot(z_axis, plane_normal));
        // 中心投影Jacobian行列式(面积比：投影后/成像平面)
        energy *= gnomonicProjection(pos, focus_center, helio_pos, d_x, d_y);
        // 以定日镜中心与像素中心连线定光线方向
        //    energy *= project_to_recv(pos, helio_pos, z_axis, plane_normal);
        atomicAdd(&(d_array_[my_id]), energy);
    }

    /**
     * 圆柱形接收器表面上定日镜光斑计算，单面定日镜，bHFLCAL
     * @param d_array_ 存储结果的数组
     * @param theta_0 圆柱表面由此角度逆时针展开
     * @param row_col d_array_数组大小
     * @param focus_center 定日镜聚焦中心
     * @param recv_size x半径y为接收器高
     * @param recv_center_pos 圆柱底面中心世界坐标
     * @param x_axis 成像平面上的x轴
     * @param y_axis 成像平面上的y轴
     * @param z_axis 成像平面上的z轴
     */
    __global__ void bHFLCALCylinderMapFlux(float* d_array_, float theta_0, float flux_param,
                                           float sigma1, float sigma2, float correlation,
                                           float3 center_bias, int2 row_col, float3 focus_center,
                                           float3 recv_size, float3 recv_center_pos, float3 x_axis,
                                           float3 y_axis, float3 z_axis)
    {
        long long my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= row_col.x * row_col.y)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        int2 pos_ij = make_int2(my_id % row_col.x, my_id / row_col.x);
        // 计算网格点的世界坐标
        float pixel_angle =
            theta_0 + 1.0F * fdividef((pos_ij.x + 0.5F), row_col.x) * 2 * inter_op::kMathPi;
        float sin_pixel_angle = 0;
        float cos_pixel_angle = 0;
        __sincosf(pixel_angle, std::addressof(sin_pixel_angle), std::addressof(cos_pixel_angle));
        float3 pos = recv_center_pos +
                     make_float3(cos_pixel_angle * recv_size.x, 0, sin_pixel_angle * recv_size.x) +
                     make_float3(0, 1.0F * fdividef((pos_ij.y + 0.5F), row_col.y) * recv_size.y, 0);
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 中心投影找对应点
        //    global_func::calcIntersection(z_axis, focus_center, pos, normalize(helio_pos-pos),
        //    projected_pos);
        // 重心偏移
        //    projected_pos -= center_bias;
        std::tie(focus_center, std::ignore) = intersect(z_axis, focus_center, center_bias, -z_axis);
        // 投影点在成像平面上的二维坐标
        // bfloat163 x_axis_t = make_bfloat163(x_axis);
        // bfloat163 y_axis_t = make_bfloat163(y_axis);
        // bfloat163 temp_pos = make_bfloat163(projected_pos - focus_center);
        // bfloat16 correlation_t = correlation;
        // bfloat16 sigma1_t = sigma1;
        // bfloat16 sigma2_t = sigma2;
        float u = dot(x_axis, projected_pos - focus_center);
        float v = dot(y_axis, projected_pos - focus_center);
        // auto u = dot(x_axis_t, temp_pos);
        // auto v = dot(y_axis_t, temp_pos);
        // 根据椭圆高斯分布计算格点处的能量
        // energy = E * f(x,y)
        // float energy = flux_param *
        //                float(hexp(bfloat16(
        //                    bfloat16(-0.5F) / (bfloat16(1) - (correlation_t * correlation_t)) *
        //                        ((u * u) / (sigma1_t * sigma1_t) + (v * v) / (sigma2_t *
        //                        sigma2_t)) -
        //                    bfloat16(2) * correlation_t * u * v / (sigma2_t * sigma1_t))));
        // float energy = flux_param *
        //                float(hexp(bfloat16(fdividef(-0.5F, (1 - (correlation * correlation))) *
        //                              (fdividef((u * u), (sigma1 * sigma1)) +
        //                               fdividef((v * v), (sigma2 * sigma2)) -
        //                               2 * correlation * u * fdividef(v, (sigma2 * sigma1))))));
        float energy = flux_param * __expf(fdividef(-0.5F, (1 - (correlation * correlation))) *
                                           (fdividef((u * u), (sigma1 * sigma1)) +
                                            fdividef((v * v), (sigma2 * sigma2)) -
                                            2 * correlation * u * fdividef(v, (sigma2 * sigma1))));
        // 乘投影系数
        // 斜平行投影Jacobian行列式
        //    energy *= parallelProjectionParam(x_axis,y_axis,z_axis,pos,helio_pos,plane_normal);
        // 中心投影Jacobian行列式(面积比：投影后/成像平面)
        // energy *= gnomonicProjection(pos,focus_center,helio_pos,d_x,d_y);
        // energy *=
        //     std::abs(dot(z_axis, make_float3(std::cos(pixel_angle), 0, std::sin(pixel_angle))));
        // //    atomicAdd(&(d_array_[myId]), energy);
        // if (dot(z_axis, make_float3(std::cos(pixel_angle), 0, std::sin(pixel_angle))) >= 0)
        // {
        //     atomicAdd(&(d_array_[my_id]), energy * fmaxf(dot(z_axis,
        //     make_float3(std::cos(pixel_angle), 0, std::sin(pixel_angle))), 0));
        // }
        // if (energy * fmaxf(dot(z_axis, make_float3(cos_pixel_angle, 0, sin_pixel_angle)), 0) > 1000)
        // {
        //     printf("u=%f v=%f sigma1=%f sigma2=%f flux_param=%f\n", u, v, sigma1, sigma2, flux_param);
        // }
        //     printf("%f\n", energy * fmaxf(dot(z_axis, make_float3(cos_pixel_angle, 0, sin_pixel_angle)), 0));

        // 斜平行投影
        atomicAdd(&(d_array_[my_id]),
                  energy * fmaxf(dot(z_axis, make_float3(cos_pixel_angle, 0, sin_pixel_angle)), 0));
    }

    __global__ void cylinderSceneSimulateV4(float3* d_image_plane_x_z_axes, float4* d_center_bias,
                                            float3* d_flux_params, int n, int length,
                                            float3 recv_size, float3 recv_pos, int2 row_col,
                                            float theta_0, float* d_array,
                                            const float* d_normal_distribution)
    {
        __shared__ float shared_normal_distribution[1024];
        shared_normal_distribution[threadIdx.x] = d_normal_distribution[threadIdx.x];
        // printf("%f\n", shared_normal_distribution[threadIdx.x]);
        __syncthreads();

        long long my_id = getThreadID();
        int helio_batch = my_id / (row_col.x * row_col.y);
        if (helio_batch * length >= n)
        {
            return;
        }
        my_id %= row_col.x * row_col.y;

        int2 pos_ij = make_int2(my_id % row_col.x, my_id / row_col.x);
        // 计算网格点的世界坐标
        float pixel_angle = theta_0 + 1.0F * (pos_ij.x + 0.5F) / row_col.x * 2 * inter_op::kMathPi;
        float3 pos =
            recv_pos +
            make_float3(cos(pixel_angle) * recv_size.x, 0, sin(pixel_angle) * recv_size.x) +
            make_float3(0, 1.0F * (pos_ij.y + 0.5F) / row_col.y * recv_size.y, 0);
        float local_total_energy = 0.0F;

        for (int helio_id = helio_batch * length;
             helio_id < n && helio_id < (helio_batch + 1) * length; helio_id++)
        {
            // // 读入定日镜信息
            float3 x_axis = d_image_plane_x_z_axes[static_cast<std::ptrdiff_t>(helio_id) * 2];
            float3 z_axis = d_image_plane_x_z_axes[static_cast<std::ptrdiff_t>(helio_id) * 2 + 1];
            float3 y_axis = cross(z_axis, x_axis);
            if (dot(z_axis, make_float3(cos(pixel_angle), 0, sin(pixel_angle))) <= 0)
            {
                continue;
            }
            // // 重心偏移
            float3 biased_center = make_float3(d_center_bias[helio_id].x, d_center_bias[helio_id].y,
                                               d_center_bias[helio_id].z);
            float sb_factor = d_center_bias[helio_id].w / 100.0f;

            // 投影点在成像平面上的二维坐标
            float sigma1 = d_flux_params[helio_id].x;
            float sigma2 = d_flux_params[helio_id].y;

            float u = abs(dot(x_axis, pos - biased_center));
            float v = abs(dot(y_axis, pos - biased_center));
            float energy = (1.0F - sb_factor) * d_flux_params[helio_id].z *
                           __expf(fdividef(-0.5F, (1 - (0 * 0))) *
                                  (fdividef((u * u), (sigma1 * sigma1)) +
                                   fdividef((v * v), (sigma2 * sigma2)) -
                                   2 * 0 * u * fdividef(v, (sigma2 * sigma1))));
            local_total_energy += energy;
            // 为了使用预计算表格，需要将N(0,a)→N(0,1)，就是在数值上除以标准差
            // float u = abs(dot(x_axis, pos - biased_center)) / sigma1;
            // float v = abs(dot(y_axis, pos - biased_center)) / sigma2;

            // int times = 100;
            // if (u < 2.5 && v < 2.5)
            // {
            //     // 做一个线性插值
            //     // 预计算表将0-4离散为0，0.01，0.02，..., 4.00
            //     // 例如u = 0.015时，它比较接近的取值应该是0.01处的预计算值和0.02处的均值
            //     auto line_alpha_u = abs(computeAlpha(u, times));
            //     auto line_alpha_v = abs(computeAlpha(v, times));
            //     // printf("%f, %f\n", line_alpha_u, line_alpha_v);

            //     auto d_normal_u =
            //         shared_normal_distribution[min(abs(int(u * times)), 1023)] *
            //             (1 - line_alpha_u) +
            //         line_alpha_u * shared_normal_distribution[min(abs(int(u * times)), 1023) +
            //         1];
            //     auto d_normal_v =
            //         shared_normal_distribution[min(abs(int(v * times)), 1023)] *
            //             (1 - line_alpha_v) +
            //         line_alpha_v * shared_normal_distribution[min(abs(int(v * times)), 1023) +
            //         1];
            //     // printf("%f, %f\n", d_normal_u, d_normal_v);

            //     auto energy = (1.0f - sb_factor) * d_flux_params[helio_id].z * 2 *
            //                   inter_op::kMathPi * d_normal_u * d_normal_v;
            //     // printf("%f\n", energy);

            //     //     // 这一步是斜平行投影
            //     energy *= std::abs(
            //         dot(z_axis, make_float3(std::cos(pixel_angle), 0, std::sin(pixel_angle))));
            //     // printf("%f\n", energy);
            //     local_total_energy += energy;
            // }
        }
        atomicAdd(&(d_array[my_id]), local_total_energy);
    }
} // namespace

SimplifiedMethod::SimplifiedMethod(float3* helio_vertices, float3* recv_vertices,
                                   float3 helio_normal, float3 recv_normal, float3 helio_pos,
                                   float3 target_center, float3 helio_size, int2 row_col, float DNI,
                                   cudaStream_t stream)
{
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        cudaStreamCreate(&stream_);
    }
    setHelioVertices(helio_vertices);
    setRecvVertices(recv_vertices);
    helio_normal_ = helio_normal;
    plane_normal_ = recv_normal;
    helio_pos_ = helio_pos;
    focus_center_ = target_center;
    helio_size_ = helio_size;
    row_col_ = row_col;
    DNI_ = DNI;
    shadow_block_factor_ = 0;
    helio_reflectivity_ = 0.88;
    center_bias_ = ::make_float3(0, 0, 0);

    CHECK_CUDA_ERRORS(
        cudaMallocAsync(std::addressof(d_array_), sizeof(float) * row_col.x * row_col.y, stream_));
    CHECK_CUDA_ERRORS(cudaMemsetAsync(d_array_, 0, sizeof(float) * row_col.x * row_col.y, stream_));
}

SimplifiedMethod::SimplifiedMethod(float3* helio_vertices, float3 recv_size, float3 helio_normal,
                                   float3 recv_pos, float3 helio_pos, float3 target_center,
                                   float3 helio_size, int2 row_col, float DNI, float* d_array,
                                   cudaStream_t stream)
{
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        cudaStreamCreate(&stream_);
    }

    setHelioVertices(helio_vertices);
    recv_size_ = recv_size;
    helio_normal_ = helio_normal;
    recv_pos_ = recv_pos;
    helio_pos_ = helio_pos;
    focus_center_ = target_center;
    helio_size_ = helio_size;
    row_col_ = row_col;
    DNI_ = DNI;
    shadow_block_factor_ = 0;
    helio_reflectivity_ = 0.88;
    center_bias_ = ::make_float3(0, 0, 0);

    if (d_array != nullptr)
    {
        this->d_array_ = d_array;
        this->external_d_array_ = true;
    }
    else
    {
        CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_array_),
                                          sizeof(float) * row_col.x * row_col.y, stream_));
        CHECK_CUDA_ERRORS(
            cudaMemsetAsync(d_array_, 0, sizeof(float) * row_col.x * row_col.y, stream_));
    }
}

static auto primeFactorization(int num) -> std::vector<int>
{
    const auto sqrt_num = std::sqrt(num);
    std::vector<int> result;
    for (int i = 2; i <= sqrt_num; i++) // 质数可以从2开始
    {
        // 判断是否为因数&是否为质数
        if (num % i == 0)
        {
            while (num % i == 0)
            {
                num /= i;
                result.emplace_back(i);
            }
        }
    }
    if (num != 1)
    {
        result.emplace_back(num);
    }
    return result;
}

SimplifiedMethod::SimplifiedMethod(float3* d_image_plane_x_z_axes, float4* d_center_bias,
                                   float3* d_flux_params, int n, float3 recv_size, float3 recv_pos,
                                   int2 row_col, cudaStream_t stream)
{
    if (stream == nullptr)
    {
        cudaStreamCreate(&stream_);
    }
    else
    {
        this->stream_ = stream;
    }
    // n = 100;

    /*1. 将数组们copy至gpu memory*/

    // 读入预计算表格
    if (d_normal_distribution == nullptr)
    {
        loadNormalDistributionArray();
    }

    /*2. 将接收器坐标转为接收器底部中心*/
    float3 cylinder_bottom_center =
        ::make_float3(recv_pos.x, recv_pos.y - 0.5F * recv_size.y, recv_pos.z);

    // 在cuda中预备存储计算结果的空间
    CHECK_CUDA_ERRORS(
        cudaMallocAsync(std::addressof(d_array_), sizeof(float) * row_col.x * row_col.y, stream_));
    CHECK_CUDA_ERRORS(cudaMemsetAsync(d_array_, 0, sizeof(float) * row_col.x * row_col.y, stream_));

    /*3. 设置线程数*/
    // 设置线程数
    const int n_threads = 1024;
    const int length = 512;
    const int batch = (n + length - 1) / length;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, row_col.x * row_col.y * batch, true);

    // std::printf("Calculation: %d\n", (long long)row_col.x * row_col.y * n);

    cylinderSceneSimulateV4<<<n_blocks, n_threads, 0, stream_>>>(
        d_image_plane_x_z_axes, d_center_bias, d_flux_params, n, length, recv_size,
        cylinder_bottom_center, row_col, 0, d_array_, d_normal_distribution);
    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
    // CHECK_CUDA_ERRORS(cudaFreeAsync(d_helio_x_z_axes, stream_));
    // CHECK_CUDA_ERRORS(cudaFreeAsync(d_flux_params, stream_));
}

SimplifiedMethod::~SimplifiedMethod()
{
    if (d_array_ != nullptr && !this->external_d_array_)
    {
        CHECK_CUDA_ERRORS(cudaFreeAsync(d_array_, stream_));
        d_array_ = nullptr;
    }
    if (d_normal_distribution != nullptr)
    {
        CHECK_CUDA_ERRORS(cudaFreeAsync(d_normal_distribution, stream_));
        d_normal_distribution = nullptr;
    }
    if (stream_ != nullptr && !this->external_stream_)
    {
        cudaStreamDestroy(stream_);
    }
}

void SimplifiedMethod::setHelioVertices(float3* helio_vertices)
{
    helio_pos_ = ::make_float3(0, 0, 0);
    std::copy(helio_vertices, helio_vertices + 4, helio_vertices_);
    for (int i = 0; i < 4; i++)
    {
        // helio_vertices_[i] = helio_vertices[i];
        helio_pos_ += helio_vertices[i];
    }
    helio_pos_ /= 4;
}

void SimplifiedMethod::setRecvVertices(float3* recv_vertices)
{
    std::copy(recv_vertices, recv_vertices + 4, recv_vertices_);
    // for (int i = 0; i < 4; i++)
    // {
    //     recv_vertices_[i] = recv_vertices[i];
    // }
}

void SimplifiedMethod::drawFluxMap()
{
    // 解析模型中定日镜镜面总能量计算
    float distance = length(focus_center_ - helio_pos_);
    float m_aa = (distance <= 1000)
                     ? (0.99321 - 0.0001176 * distance + 1.97 * 1e-8 * distance * distance)
                     : std::exp(-0.0001106 * distance);
    float s = helio_size_.x * helio_size_.z;
    float cos_w = dot(normalize(focus_center_ - helio_pos_), helio_normal_);

    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[2];
    for (int i = 1; i < 3; i++)
    {
        std::tie(interv[i - 1], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[1] - interv[0]);
    y_axis = cross(z_axis, x_axis);

    float flux_param = peak_ * (1.0F - shadow_block_factor_);

    // d_x, d_y
    float3 d_x = (recv_vertices_[1] - recv_vertices_[0]) / row_col_.y;
    float3 d_y = (recv_vertices_[3] - recv_vertices_[0]) / row_col_.x;
    float3 recv_vertex_0 = recv_vertices_[0] + d_x / 2 + d_y / 2;
    //    float3 recv_vertex[] = {make_float3(10,170,0), make_float3(10,190,0),
    //    make_float3(-10,170,0)}; float3 d_x = (recv_vertex[1]-recv_vertex[0]) / row_col_.y; float3
    //    d_y = (recv_vertex[2]-recv_vertex[0]) / row_col_.x; float3 recv_vertex_0 = recv_vertex[0]
    //    + d_x / 2 + d_y / 2;

    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, row_col_.x * row_col_.y, true);
    //
    bHFLCALMapFlux<<<n_blocks, n_threads, 0, stream_>>>(
        d_array_, flux_param, sigma1_, sigma2_, correlation_, row_col_, focus_center_, center_bias_,
        recv_vertex_0, d_x, d_y, x_axis, y_axis, z_axis, plane_normal_, helio_pos_);

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

void SimplifiedMethod::drawFluxMapCylinder()
{
    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[2];
    for (int i = 1; i < 3; i++)
    {
        std::tie(interv[i - 1], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[1] - interv[0]);
    y_axis = cross(z_axis, x_axis);

    float flux_param = peak_ * (1.0F - shadow_block_factor_);

    // solar::println("Before kernel bHFLCALCylinderMapFlux");
    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, row_col_.x * row_col_.y, true);
    // solar::println(solar::stringConcatenate("n_blocks: ", std::to_string(n_blocks.x), ", ",
    //                                         std::to_string(n_blocks.y), ", ",
    //                                         std::to_string(n_blocks.z)));
    //
    float3 cylinder_bottom_center =
        ::make_float3(recv_pos_.x, recv_pos_.y - 0.5F * recv_size_.y, recv_pos_.z);
    bHFLCALCylinderMapFlux<<<n_blocks, n_threads, 0, stream_>>>(
        d_array_, 0, flux_param, sigma1_, sigma2_, correlation_, center_bias_, row_col_,
        focus_center_, recv_size_, cylinder_bottom_center, x_axis, y_axis, z_axis);

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

auto SimplifiedMethod::loadNormalDistributionArray() -> bool
{
    try
    {
        const auto* array_ptr = &normal_distribution_array_precomp[0];
        CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(this->d_normal_distribution),
                                          sizeof(float) * 1024, stream_));
        CHECK_CUDA_ERRORS(
            cudaMemsetAsync(this->d_normal_distribution, 0, sizeof(float) * 1024, stream_));
        cpu2gpuAsync(this->d_normal_distribution, array_ptr,
                     std::size(normal_distribution_array_precomp), stream_);
        return true;
    }
    catch (std::exception& e)
    {
        throw e;
    }
}

void SimplifiedMethod::setCudaStream(const cudaStream_t& stream)
{
    // if (this->stream_ != nullptr && !this->external_stream_)
    // {
    //     cudaStreamDestroy(this->stream_);
    // }
    this->stream_ = stream;
}
